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Point-based  surface  processing  has  developed  into  an  attractive  alternative  to  mesh-based  process¬ 
ing  techniques  for  a  number  of  geometric  modeling  applications.  By  working  with  point  clouds 
directly,  any  processing  is  based  on  the  given  raw  data  and  its  underlying  geometry  rather  than 
any  arbitrary  intermediate  representations  and  generally  artificial  connectivity  relations.  We  ex¬ 
tend  this  principle  into  the  area  of  subdivision  surfaces  by  introducing  the  notion  of  meshless,  or 
point  cloud,  geometric  subdivision.  Our  meshless  subdivision  approach  replaces  the  role  of  mesh 
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the  computation  of  intrinsic  means  on  euclidean  manifolds. 
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1.  INTRODUCTION 

Mesh  subdivision  has  developed  into  a  powerful  and  widely-used  tool  for  the  free¬ 
form  design,  editing  and  representation  of  smooth  surfaces.  Subdivision  schemes 
recursively  apply  a  local  subdivision  operator  to  a  coarse  base  mesh  thereby  produc¬ 
ing  a  sequence  of  refined  meshes  which  quickly  converges  to  a  smooth  limit  surface. 
The  advantages  of  mesh  subdivision  include  guaranteed  global  surface  smoothness 
whilst  supporting  local  feature  control,  the  ability  to  handle  surfaces  of  arbitrary 
topology  and  being  efficiently  and  simple  to  apply  once  a  base  mesh  is  available 
[Dyn  and  Levin  2002;  Zorin  and  Schroder  2000]. 

Unfortunately,  when  dealing  with  point-sampled  geometry,  mesh  subdivision  re¬ 
quires  frequently  costly  and  generally  non-geometric  surface  reconstruction,  see, 
e.g.,  [Amenta  et  al.  1998;  Amenta  et  al.  2001;  Bernadini  et  al.  1999;  Boissonnat 
and  Cazals  2000;  Curless  and  Levoy  1996;  Edelsbrunner  and  Miicke  1994],  often 
followed  by  mesh  simplification  ([Gotsman  et  al.  2002]  and  references  therein), 
parameterisation  [Floater  and  Hermann  2004]  and  remeshing  ([Alliez  et  al.  2003] 
and  references  therein),  all  pre-processing  steps  to  obtain  a  base  mesh.  During 
the  reconstruction  step,  any  measurement  noise,  misalignment  of  scans,  etc.  may 
translate  into  topological  artifacts  in  the  form  of  erroneous  connectivity  and  genus 
[Wood  et  al.  2004] .  This  hinders  subsequent  mesh  simplification  and  remeshing  and 
thus  mesh  subdivision  processing.  Also,  in  the  case  of  extremely  high-dimensional 
manifolds  by  samples  [Tenenbaum  et  al.  2000],  mesh  pre-processing  breaks  down 
at  the  surface  reconstruction  step  and  it  needs  to  be  worked  directly  with  the  point 
cloud. 

This  paper  advocates  the  use  of  meshless,  or  point  cloud,  geometric  subdivision. 
We  propose  to  avoid  the  consideration  of  mesh  connectivity  graphs  and  the  asso¬ 
ciated  pre-processing  steps  and  instead  to  work  with  the  point-sampled  geometry 
directly  using  intrinsic  subdivision  rules.  In  this  paper,  we  show  the  conceptual  via¬ 
bility  of  this  notion  of  meshless  geometric  subdivision  by  introducing  a  first  intrinsic 
meshless  subdivision  scheme  using  geodesic  centroids  of  intrinsic  point  neighbour¬ 
hoods.  We  put  forward  a  new  method  for  the  computation  of  these  geodesic  means, 
which  by  itself  is  of  interest  in  other  areas  such  as  intrinsic  statistical  shape  analysis 
[Fletcher  et  al.  2003]  and  variational  theory  [dost  1994;  Kendall  1990]. 

Some  of  the  related  geometric  operations  may  be  approximated  in  an  Euclidean 
context  when  working  with  large  sampling  densities,  regular  meshes  and  very  local 
subdivision  rules.  Working  with  the  point-sampled  geometry  directly  and  intrinsi¬ 
cally,  however,  avoids  the  need  for  non-geometric  pre-processing  steps  and  special 
rules  dealing  with  irregular  mesh  connectivity.  Unlike  the  non-geometric  nature  of 
Euclidean-based  (pre-)processing,  meshless  intrinsic  subdivision  inherently  captures 
the  non-linear  intrinsic  structure  of  the  object  geometry  and  the  principle  applies 
equally  well  to  three  and  higher-dimensional  surfaces.  Although  it  is  possible  to 
obtain  a  geodesic  mesh  representation  of  a  point-sampled  surface  [Peyre  and  Cohen 
2003;  Sifri  et  al.  2003],  to  continue  performing  subdivision  truly  intrinsically,  this 
mesh  would  have  to  be  modified  repeatedly  during  each  iteration  to  re-enforce  its 
intrinsic  nature  and  more  than  one  mesh  would  be  required  to  be  able  to  compute 
correct  geodesic  distances. 

Our  algorithm  operates  intrinsically  throughout  without  the  need  for  prior  surface 
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reconstruction  with  the  help  of  the  geodesic  distance  mapping  algorithm  for  point 
clouds  presented  in  Memoli  and  Sapiro  [2003].  Following  an  overview  of  related 
work  in  Section  2  and  some  mathematical  preliminaries  in  Section  3,  we  therefore 
summarise  the  main  aspects  of  this  technique  in  Section  4.  Section  5  presents  our 
intrinsic  meshless  subdivision  algorithm.  Section  6  gives  experimental  results  and 
comments  on  implementational  aspects.  In  the  concluding  Section  7,  we  briefly 
remark  on  our  ideas  and  some  preliminary  results  for  the  theoretical  analysis  of 
intrinsic  meshless  subdivision  schemes. 


2.  RELATED  WORK 

We  start  with  a  summary  of  the  ideas  underpinning  mesh-based  subdivision  of 
surfaces  in  The  overview  is  motivational  in  nature,  for  a  thorough  formal 
treatment  of  mesh  subdivision,  see  Dyn  and  Levin  [2002].  The  section  is  concluded 
with  remarks  on  recent  progress  in  meshless,  point-based  surface  processing  related 
to  our  work. 

Following  the  notation  of  Dyn  and  Levin  [2002],  surface  subdivision  schemes  con¬ 
sist  of  a  subdivision  operator  S  recursively  applied  to  control  nets  Ni  =  N{V\  ,  F^) 
of  arbitrary  topology,  with  I  G  Zq  denoting  the  subdivision  level,  V  representing 
a  set  of  control  vertices  in  and  E  and  F  describing  the  topological  relations  in 
the  form  of  edges  and  faces  respectively.  The  iterative  application  of  this  scheme 
generates  a  sequence  A^/+i  =  SNi.  More  speciflcally,  starting  with  a  coarse  base  net 
Nq,  at  each  iteration,  new  control  vertices  are  inserted  and  connected  according  to 
the  scheme’s  refinement  rule  and  re-positioned  following  the  operator’s  geometric 
averaging  rule.  Both  the  refinement  and  the  geometric  averaging  rule  give  the  po¬ 
sition  of  control  vertices  in  Ni^i  in  the  form  of  weighted  averages  of  topologically 
neighbouring  vertices  in  Ni .  The  careful  choice  of  these  rules  in  relation  to  the  con¬ 
trol  vertex  valency,  i.e.  the  number  of  edges  emanating  from  the  vertex,  guarantees 
the  convergence  of  the  scheme,  in  each  component  and  in  the  uniform  norm,  to  a 
limit  surface  of  provable  continuity. 

Not  every  existing  mesh  subdivision  operator  allows  for  this  simple  distinction 
between  a  topological  refinement  and  a  geometric  averaging  rule.  However,  those 
that  do  allow  for  this  kind  of  distinction,  include  the  most  widely-used  schemes. 
For  example,  the  Loop  [1987]  subdivision  scheme  for  triangular  control  nets  may 
be  cast  in  this  form.  In  the  case  of  this  scheme,  the  refinement  rule  adds  points 
related  to  each  edge  in  the  triangulation  using  face  splitting.  The  averaging  rule  of 
points,  which  depends  on  the  vertices  and  edges  of  the  non-reflned  triangulation, 
then  yields  the  final  positions  of  the  vertices  in  the  refined  triangulation. 

Point-based  surface  processing  advocates  the  use  of  points  as  editing  and  display 
primitive  thereby  avoiding  the  need  for  a  mesh  connectivity  graph  and  the  manifold 
reconstruction  step.  Recent  progress  in  this  held  has  been  sufficiently  substantial 
to  realise  powerful  point-based  editing,  free-form  and  multiresolution  modelling 
[Alexa  et  al.  2003;  Linsen  2001;  Pauly  et  al.  2002;  Pauly  et  al.  2003;  Zwicker 
et  al.  2002]  and  visualisation  [Grossman  and  Dally  1998;  Kalaiah  and  Varshney 
2001;  Levoy  and  Whitted  1985;  Pfister  et  al.  2000;  Rusinkiewicz  and  Levoy  2000; 
Zwicker  et  al.  2001;  2002]  alternatives  to  mesh-based  processing  methods.  With 
numerous  applications  in  medical  imaging,  reverse  engineering,  cultural  heritage 
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preservation,  natural  science  experimentation,  geoscience,  numerical  simulations, 
etc.  acquiring  point  sets  of  significant  density  from  three  or  higher-dimensional 
manifold  surfaces,  it  is  particularly  attractive  to  avoid  the  generally  costly,  error- 
prone  and  non-geometric  mesh-related  pre-processing  steps  and  to  work  with  point- 
sampled  geometry  directly. 

In  this  paper,  we  propose  to  replace  the  role  of  mesh  connectivity  in  subdivision 
with  intrinsic  point  neighbourhood  information  and  formulate  a  set  of  meshless 
refinement  and  geometric  averaging  rules  in  the  form  of  weighted  geodesic  centroids 
of  these  local  neighbourhoods. 

To  the  best  of  our  knowledge,  only  Fleishman  et  al.  [2003]  and  Guennebaud  et  al. 
[2004]  touch  upon  the  notion  of  point  cloud  subdivision.  In  Fleishman  et  al.  [2003], 
the  authors  generate  progressive  levels-of-detail  of  point  clouds  by  transferring  the 
mesh-based  idea  of  subdivision  displacement  maps  to  the  point  cloud  case.  They 
devise  a  point  cloud  simplification  method  for  the  generation  of  a  base  point  set  and 
present  both  a  projection  and  a  local,  uniform  upsampling  operator  with  the  help 
of  local  surface  reconstruction  using  Moving  Least  Squares  [Alexa  et  al.  2003].  The 
authors  successfully  mimic  the  principle  of  mesh-based  subdivision  displacement 
mapping  for  surfaces  in  point  cloud  form  but  do  not  consider  the  idea  of  meshless 
subdivision. 

Similarly,  Guennebaud  et  al.  [2004]  ^  are  concerned  with  the  upsampling  of  “surfel 
sets”,  i.e.  input  points  equipped  with  normal,  local  sampling  density  (surfel  radius) 
and  texture  information,  for  magnified  point  rendering  by  splatting.  Their  inter- 
polatory  method  requires  that  the  underlying  point  cloud  is  regularly  uniformly 
distributed  and  noise-free  and  that  features  such  as  crease  lines  have  been  detected 
and  adequately  sampled  in  a  pre-processing  step.  Their  technique  is  restricted  in 
applicability  to  surfaces  in  and  intended  for  the  operation  on  top  of  a  splat- 
ting  algorithm  such  as  [Zwicker  et  al.  2001;  2002]  providing  the  surfel  information. 
As  is  typically  the  case  for  independently  determined  neighbourhoods  such  as  the 
authors’  polygonal  fan  neighbourhoods,  their  (extrinsic)  proximity  and  refinement 
operators  are  affected  by  overlapping  neighbour  relations  and  the  need  for  refine¬ 
ment  rules  varying  with  the  number  of  neighbours  of  the  point  under  consideration. 
This  interpolatory  mesh  subdivision-inspired  method  is  used  for  the  locally  smooth 
upsampling  of  surfel  sets  as  part  of  a  dynamic  splatting  algorithm,  the  more  general 
notion  of  meshless  subdivision  is  not  considered. 

This  paper  is  concerned  with  the  notion  of  meshless  surface  subdivision  and 
makes  no  assumptions  on  the  availability  of  normal  or  local  density  information  or 
the  regularity  of  the  input  point  cloud.  An  intrinsic  neighbourhood  concept  is  used 
which  is  based  on  a  partitioning  of  the  surface  and  avoids  the  shortcomings  of  point 
neighbourhoods  determined  independently  from  each  other.  The  approach  is  not 
restricted  to  surfaces  in  but  extends  to  higher  dimensions.  Inspired  by  recent 
work  on  geodesic  curve  subdivision  [Wallner  and  Dyn  2003;  2004;  2005],  it  offers  a 
theoretical  basis  for  global  convergence  and  smoothness  analysis. 

In  the  following,  we  present  our  intrinsic  framework  for  meshless  subdivision. 
We  start  with  notation  and  the  definition  of  some  key  notions  used  throughout  the 


^This  work  was  produced  at  about  the  same  time  as  ours,  see  [Moenning  et  al.  2004]  for  an 
abbreviated  publication  of  the  ideas  we  now  put  forward  in  detail. 
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paper.  This  is  followed  by  the  discussion  of  the  intrinsic  distance  mapping  algorithm 
of  Memoli  and  Sapiro  [2003]  and  the  presentation  of  the  framework  itself. 

3.  PRELIMINARIES 

Let  M  represent  a  differentiable  (smooth),  compact  and  connected  Riemannian 

submanifold  in  K™,  m  >  3.  By  intrinsic  processing,  we  mean  processing  directly 

on  this  submanifold  rather  than  in  its  embedding  space.  M  is  represented  by  a 

finite  point  set,  or  (unstructured)  point  cloud,  P  =  {pi,P2,  ■  ■  ■  :Pn}  C  M.  The 

Riemannian  metric  on  M  at  point  x  G  M  is  a,  smoothly  varying  inner  product  (., .) 

1 

on  the  tangent  space  T^M.  The  norm  of  a  vector  v  in  TxM  is  given  by  ||u||  =  {v,  v)  ^ . 
We  endow  M  with  the  metric  inherited  from  M™,  hence  {v,w)  will  coincide  with 
the  usual  inner  product  for  vectors  v  and  w  in  M™.  Consider  a  (sectionally)  smooth 
curve  7  :  [a,  5]  G  M  ^  M  parameterised  by  t.  The  length  of  7(t)  follows  from 
integrating  the  norm  of  its  tangent  vectors,  j{t),  along  the  curve,  i.e. 

Hl)=[  \\i{t)\\dt 

J  a 

The  curve  7  is  called  the  (minimising)  geodesic  from  a  point  x  to  a  point  y  in  M,  if 
it  represents  the  minimum  length-curve  among  all  the  curves  on  M  joining  x  and 
y.  Since  we  are  assuming  M  to  be  compact,  it  is  geodesically  complete  and  there 
exists  at  least  one  such  curve  on  M .  This  curve  may  not  be  unique.  The  length 
of  the  geodesic  between  x  and  y  gives  the  intrinsic,  or  geodesic,  distance,  dM{x,y), 
between  the  points,  i.e. 

dM{x,y)  =  inf{L(7)}, 

7 

with  7(0)  =  X  and  7(6)  =  y.  The  function  giving  the  intrinsic  distance  from  a 
point  X  G  M  to  every  point  in  M,  dM{x, .),  is  called  the  intrinsic  distance  function, 
or  intrinsic  distance  map,  of  x.  Unique  geodesics  between  two  points  x  and  y  on 
M  may  thus  be  computed  from  dM{x,  •)  by  backtracking  from  y  towards  x  in  the 
direction  of  the  (negative)  gradient  of  dM{x,-). 

The  extrinsic  distance  between  points  x  and  y  on  M,  d{x,y),  is  computed  in  the 
metric  of  the  embedding  space.  Since  we  are  concerned  with  manifold  surfaces  in 
M™,  the  extrinsic  distance  is  Euclidean  and  given  by  the  length  of  the  Euclidean 
line  segment  between  x  and  y  in  the  ambient  space.  Apart  from  its  end  points, 
this  line  segment  generally  does  not  lie  on  the  manifold.  More  detail  on  the  above 
notions  may  be  found  in,  for  example,  Chavel  [1997]. 

By  intrinisic,  or  geodesie,  centroid,  we  understand  the  mean  of  a  local  neighbour¬ 
hood  of  points  on  M  computed  in  terms  of  intrinsic  distances  between  the  points. 
This  is  to  be  distinguished  from  the  extrinsic  centroid  of  the  subset,  computed  us¬ 
ing  Euclidean  distances  in  the  ambient  space  and  subsequently  projected  onto  the 
manifold  surface. 

Our  meshless  subdivision  scheme  makes  extensive  use  of  the  geodesic  Voronoi 
diagram  concept:  Define  the  bisector  BS{pi,Pj)  of  Pi,  Pj  G  P,  Pi  ^  pj,  as  geodesi¬ 
cally  equidistantial  loci  with  respect  to  Pi,Pj,  i.e.  BS{pi,Pj)  =  {q  G  M\dM{Pi,  q)  = 
dniPj,  ?)}•  Let  the  dominance  region  of  pi,  D{pi,pj),  denote  the  region  of  M  con¬ 
taining  Pi  bounded  by  BS{pi,pj).  The  Voronoi  region  of  pi  with  respect  to  point 
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set  P  is  given  by  R{pi,P)  =  Clp.^p p.^p.  D{pi,pj)  and  consists  of  all  points  for 
which  the  geodesic  distance  to  pi  is  smaller  than  the  geodesic  distance  to  any  other 
Pj  G  P.  The  boundary  shared  by  a  pair  of  Voronoi  regions  is  called  a  Voronoi 
edge.  Voronoi  edges  meet  at  Voronoi  vertices.  The  geodesic  Voronoi  diagram  of  P 
is  defined  as 


VD{P)=  U  dR{p,,P), 

Pi&p 

where  dR{pi,P)  denotes  the  boundary  of  R{pi,P).  Kunze  et  al.  [1997]  and  Leibon 
and  Letscher  [2000]  give  further  details  on  the  notion  of  geodesic  Voronoi  diagrams. 

4.  INTRINSIC  DISTANCE  MAPPING  ACROSS  POINT  CLOUDS 

To  compute  intrinsic  distance  maps,  and  from  them  (minimal)  geodesics,  across 
point  clouds,  we  use  the  extension  of  the  original  Fast  Marching  method  [Helmsen 
et  al.  1996;  Sethian  1996;  Tsitsiklis  1995]  to  surfaces  in  point  cloud  form  introduced 
by  Memoli  and  Sapiro  [2003].  In  the  following,  we  summarise  the  main  aspects  of 
this  technique.  For  full  details,  see  Memoli  and  Sapiro  [2001;  2003]. 

Fast  Marching  generally  represents  a  very  efficient  technique  for  the  solution  of 
front  propagation  problems  which  can  be  formulated  as  boundary  value  partial 
differential  equations.  Take  the  simple  case  of  a  front  propagating  across  a  3D  Eu¬ 
clidean  domain  with  speed,  or  weight,  function  F{v)  away  from  a  source  (boundary) 
u,  with  u  and  v  representing  3-tuples  in  We  are  interested  here  in  the  arrival 
time,  or  offset  distance,  d{u,v),  of  the  front  at  grid  point  v.  The  relationship  be¬ 
tween  the  magnitude  of  the  distance  map’s  gradient,  V(i(u,  v),  and  the  given  weight 
P(v)  at  V  can  be  expressed  as  the  following  boundary  value  formulation 

|Vd(M,u)l  =  F(v), 

with  boundary  condition  d(u,u)  =  0.  The  problem  of  determining  a  Euclidean 
weighted  distance  map  has  therefore  been  transformed  into  the  problem  of  solv¬ 
ing  a  particular  type  of  static  Hamilton- Jacobi  partial  differential  equation,  the 
nonlinear  Eikonal  equation.  For  F{v)  >  0,  this  type  of  equation  can  be  solved  ef¬ 
ficiently,  in  computationally  optimal  time,  using  conventional  Cartesian  numerics, 
see,  [Helmsen  et  al.  1996;  Sethian  1996;  Tsitsiklis  1995]. 

Memoli  and  Sapiro  [2003]  extend  the  applicability  of  this  Fast  Marching  idea  to 
the  case  of  general  co-dimension  submanifolds  in  point  cloud  form  in  three  or  higher 
dimensions.  For  simplicity,  consider  the  constant  radius  r-offset  flp,  i.e.  the  union 
of  Euclidean  balls  with  radius  r  centred  at  points  pi  &  P  (Figure  1) 

n 

rip  :=  (J  B{pi,r)  =  {cc  G  :  d{pt,x)  <  r} 

i=l 

To  approximate  the  weighted  intrinsic  distance  map,  dM{q,‘),  originating  from 
a  source  point  q  G  M  on  M,  Memoli  and  Sapiro  [2003]  suggest  computing  the 
Euclidean  distance  map  in  f2p,  denoted  dnj.-  That  is 

|VMdM(p,-)l  =  F{p), 
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Fig.  1.  Intrinsic  distance  mapping  using  Fast  Marching  for  point  clouds  operates  in  an  (not 
necessarily  constant  radius)  offset  band  consisting  of  the  union  of  balls  B{pi,  r)  centred  at  (black) 
points  Pi  of  the  surface  M  (left).  Only  those  (blue)  grid  points  falling  inside  the  offset  band 
are  considered  during  processing.  A  cross-sectional  view  of  a  constant  radius  offset  band  for  the 
Michelangelo  Youthful  data  set  is  shown  on  the  right. 

for  p  S  M  and  with  boundary  condition  dM{q,  (?)  =  0  is  approximated  by 

|Vdo^,(p,-)l  =  Art,  (2) 

for  p  G  QJp  and  boundary  condition  =  Q.  F  represents  the  (smooth) 

extension  of  the  propagation  speed  F  on  M  into  fip;  Vm  denotes  the  intrinsic 
gradient  operator.  The  problem  of  computing  an  intrinsic  weighted  distance  map 
is  therefore  transformed  into  the  problem  of  computing  an  Euclidean,  or  extrinsic, 
weighted  distance  map  in  the  offset  band  flp  around  the  surface,  i.e.  in  an  Euclidean 
manifold  with  boundary. 

Memoli  and  Sapiro  [2003]  prove  uniform  (probabilistic)  convergence  between 
these  two  distance  maps,  and  geodesics  computed  from  them,  for  both  noise-free 
and  noisy  (provided  noise  is  bounded  from  above  by  r),  randomly-sampled  point 
clouds  and  thus  show  that  the  approximation  error  between  the  intrinsic  and  ex¬ 
trinsic  distance  maps  is  of  the  same  theoretical  order  as  that  of  the  conventional 
Fast  Marching  algorithm.  Fast  Marching  can  therefore  be  used  to  approximate  the 
solution  to  (2)  in  a  computationally  optimal  manner  and  without  the  need  for  any 
prior  surface  reconstruction  by  only  slightly  modifying  the  conventional  Cartesian 
Fast  Marching  technique  to  deal  with  bounded  spaces.  This  is  achieved  by  simply 
restricting  the  grid  points  visited  by  the  conventional  Fast  Marching  algorithm  to 
those  located  in  flp.  By  performing  the  computations  within  this  offset  band,  this 
method  is  relatively  robust  in  the  presence  of  noisy  point  samples,  especially  when 
compared  to  graph-based  distance  mapping  algorithms  such  as  [Giesen  and  Wagner 
2003;  Tenenbaum  et  al.  2000]  in  which  case  the  geodesics  pass  through  the  noisy 
samples  rather  than  an  union  of  Euclidean  balls  centred  at  the  input  points. 

The  complexity  of  this  algorithm  is  0{N  log  N),  where  N  represents  the  number 
of  grid  points  located  in  the  (narrow)  offset  band  flp.  Memory  efficiency  is  achieved 
by  storing  these  grid  points  only.  Note  in  this  context  that  subject  to  the  bounds 
given  in  Memoli  and  Sapiro  [2003] ,  r  will  generally  be  small  and  does  not  have  to  be 
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constant  but  may  vary  with  each  pi.  Implementational  details  are  given  in  Section 

6.2. 

This  Fast  Marching  technique  underpins  our  geodesic  computation  algorithm 
presented  as  part  of  our  intrinsic  meshless  subdivision  framework  in  the  following 
section. 

5.  INTRINSIC  MESHLESS  SUBDIVISION 

Subdivision  schemes  incorporate  refinement  and  geometric  averaging  rules  in  the 
form  of  weighted  averages  of  local  neighbourhoods.  Mesh-based  subdivision  schemes 
are  derived  in  a  parametric  setting  ignoring  the  geometric  embedding  of  the  points 
in  space.  As  a  result,  they  are  formulated  in  terms  of  local  mesh  connectivity  rather 
than  object  geometry.  We  determine  local  neighbourhood  relations  from  intrinsic 
point  proximity  information  instead. 

We  start  this  section  with  the  discussion  of  a  suitable  intrinsic  neghbourhood 
concept.  This  is  followed  by  the  presentation  of  our  meshless  subdivision  scheme 
and  a  new  method  for  the  computation  of  geodesic  centroids,  which  is  at  the  heart 
of  this  scheme. 

5.1  Intrinsic  point  proximity  information 

Depending  on  the  acquisition  technique  used,  point-sampled  geometry  may  feature 
(locally)  non-uniform  point  distributions.  In  this  setting,  naive  neighbourhood 
concepts  such  as  simple  ball  or  k  nearest  neighbourhoods  tend  to  be  skewed,  i.e. 
the  neighbours  qj  G  Pi  are  frequently  no  longer  distributed  spherically,  i.e.  all 
around,  the  point  p  €  Pi  under  consideration  [Floater  and  Reimers  2001;  Linsen 
2001]. 

To  allow  for  local  sampling  non-uniformities,  Linsen’s  [2001]  enhanced  k  nearest 
neighbourhood,  which  enforces  a  maximum  angle  between  successive  neighbours 
around  p,  or  Floater  and  Reimer’s  [2001]  local  Delaunay  neighbourhood  may  be 
used.  Their  algorithm  projects  a  local  ball  neighbourhood  of  p  and  p  itself  into 
their  least  squares  plane  and  computes  the  planar  Delaunay  triangulation  of  the 
projected  points.  The  points’  neighbour  relations  in  this  triangulation  are  then 
taken  as  the  neighbour  relations  of  p  on  the  manifold.  We  suggest  to  collect  local 
proximity  information  by  considering  the  set  of  (Voronoi)  neighbours,  Np,  of  p  in 
the  geodesic  Voronoi  diagram  of  Pi,  VD{Pi),  instead,  i.e.  an  intrinsic  “natural 
neighbourhood”  [Sibson  1980], 

Np  =  {q  ■.  p  and  q  are  neighbours  in  VD{Pi)}, 

for  p,  q  &  Pi,  p  q  (Figure  2).  Due  to  its  intrinsic  nature  and  unlike  (Euclidean) 
ball  neighbourhoods,  points  from  disjoint  parts  of  a  surface  are  prevented  from  be¬ 
ing  assigned  to  the  same  natural  neighbourhood.  The  need  for  determining  local 
ball  radii  and  the  problem  of  locally  ill-defined  least  squares  fits  in  the  case  of  k 
nearest  neighbourhoods  do  not  arise.  Since  intrinsic  natural  neighbourhoods  are 
derived  from  a  partitioning  of  the  surface,  the  problem  of  overlapping  neighbour 
relations  is  avoided  as  well.  This  problem  is  typically  encountered  when  using 
point  neighbourhoods  computed  independently  from  each  other  such  as  k  nearest 
neighbours  and  its  variations  [Guennebaud  et  al.  2004;  Linsen  2001]  and  occurs 
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Fig.  2.  As  illustrated  here  for  the  planar  case,  the  intrinsic  Voronoi  diagrams  (solid  lines)  of 
various  moderately  irregularly  distributed  point  sets  provide  point  proximity  relations  in  the  form 
of  intrinsic  natural  neighbourhoods  (blue  points)  spherically  distributed  around  the  (red)  point 
under  consideration. 

irrespective  of  the  regularity  of  the  point  set  distribution.  Intrinsic  natural  neigh¬ 
bourhoods  are,  however,  less  straightforward  to  compute.  We  discuss  our  approach 
in  detail  in  Section  6. 

We  use  intrinsic  natural  neighbourhood  information  for  the  computation  of  local 
geodesic  centroids  as  part  of  our  meshless  subdivision  scheme  presented  next. 

5.2  An  intrinsic  meshless  subdivision  scheme 

Within  our  meshless  geometric  subdivision  framework,  geodesic  centroids  of  the 
intrinsic  neighbourhoods  discussed  in  the  previous  section  are  used  to  define  mesh¬ 
less  subdivision  rules.  More  specifically,  we  propose  the  following  set  of  rules  to  be 
applied  at  each  iteration: 


Refinement  rule:  For  each  neighbour  G  Np^  consider  the  union  of  intrinsic 
neighbours  of  p,  qj  G  Pi ,  Npq. .  Upsample  Pi  by  inserting  the  weighted  geodesic 
centroid,  c{Npq.)  G  Pi+i,  of  Npq.. 

Geometric  averaging  rule:  Replace  p  G  P;  by  the  weighted  geodesic  centroid, 
c{Np)  G  Pi+i,  of  its  intrinsic  Voronoi  neighbourhood  Np. 


This  use  of  weighted  centroids  in  the  refinement  and  geometric  averaging  rules  is 
reminiscent  of  both  classical  subdivision  schemes  [Zorin  and  Schroder  2000]  (Section 
2)  and  the  “repeated  averaging”  approach  towards  the  generation  of  subdivision 
surfaces  ([Lane  and  Riesenfeld  1980;  Oswald  and  Schroder  2003]  and  references 
therein).  These  subdivision  rules  are  incorporated  into  our  meshless  subdivision 
algorithm  as  summarised  in  Algorithm  1. 

In  its  initialisation  phase,  the  algorithm  buckets  the  input  point  set  P;  in  a  Carte¬ 
sian  grid,  subsequently  used  to  support  weighted  distance  mapping  and  geodesic 
centroid  computation.  Initialisation  is  completed  with  the  computation  of  the  (dis¬ 
crete)  geodesic  Voronoi  diagram  of  Pi,  VD{Pi).  Main  processing  loops  over  all 
points  in  Pi  and  proceeds  with  the  upsampling  of  P;  to  Pj+i  using  joint  neighbour¬ 
hood  information  readily  available  from  VD(Pi)  and  our  refinement  rule.  Following 
this  refinement  step,  VD{Pi+i)  is  computed.  Meshless  subdivision  is  concluded  with 
the  application  of  our  geometric  averaging  rule  to  each  point  in  the  refined  point 
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Input:  Point  cloud  Pi  £  R*". 

Output:  Subdivided  point  cloud  Pi+i. 

0  ***  INITIALISATION  *** 

1  Bucket  the  base  point  cloud  Pi  in  a  m-dimensional  Cartesian  grid; 

2  Compute  the  discrete  geodesic  Voronoi  diagram,  VD{Pi),  of  Pi. 

3 

4  ***  MAIN  PROCESSING  *** 

5  FOR  each  point  pi  £  Pr, 

6  Determine  the  Voronoi  neighbourhood  Np^  from  VD{Pi); 

7  FOR  each  neighbour  qj  £  Np. ; 

8  Determine  the  joint  Voronoi  neighbourhood  from  VD{Pi); 

9  Compute  the  weighted  geodesic  centroid  c{Np^q.)-, 

10  {Refinement  rule)  Upsample  Pi  to  Pi+i  by  inserting  c{Np.qj); 

11  ENDFOR 

12  ENDFOR 

13  Compute  l/D(P;+i); 

14  FOR  each  point  pi  £  Pi+i ; 

15  Determine  the  Voronoi  neighbourhood  Np.  from  UD(p+i); 

16  Compute  the  weighted  geodesic  centroid  c{Np.); 

17  {Geometric  averaging  rule)  Replace  pi  in  Pi+i  with  c{Npfi)-, 

18  ENDFOR 

Alg.  1:  One  iteration  of  meshless  subdivision  in  pseudocode. 


set  Pi+i  yielding  the  final,  subdivided  point  set. 

The  applicability  of  this  subdivision  algorithm  does  not  depend  on  a  preceding 
simplification  step.  Potential  algorithm  applications,  however,  include  the  case 
in  which  it  may  be  desirable  to  simplify  an  excessively  dense  point  cloud  P  to  a 
suitable  base  point  set  Pq  in  the  expectation  that  recursive  subdivision  of  Pq  results 
in  a  smoother,  more  regular  and  more  compact  approximation  of  the  underlying 
surface  than  given  by  P.  Since  our  subdivision  algorithm  requires  the  computation 
of  geodesic  centroids  across  the  base  point  set  Pq,  for  these  centroids  to  be  well- 
defined,  any  simplification  of  P  needs  to  be  performed  subject  to  a  minimum  point 
density  in  Pq.  Due  to  the  method’s  simple  control  of  a  guaranteed  point  density, 
its  purely  intrinsic  operation,  its  close  relationship  with  the  natural  neighbourhood 
concept  of  Section  5.1  and  its  efficient  implementability  using  the  geodesic  distance 
mapping  method  of  Memoli  and  Sapiro  [2003]  (Section  4),  we  use  the  algorithm  of 
Moenning  and  Dodgson  [2004]  for  the  simplification  of  an  input  point  cloud  P  to 
a  base  point  set  Pq  still  sufficiently  dense  to  support  the  computation  of  geodesic 
centroids.  As  another  result  of  this  pre-processing  step,  the  (discrete)  geodesic 
Voronoi  diagram  of  Pq  becomes  available  [Moenning  and  Dodgson  2004]  so  that  it 
does  not  need  to  be  computed  in  our  subdivision  algorithm’s  initialisation  phase 
and  the  natural  neighbours  of  a  point  pi  £  Pq  are  readily  known. 

By  performing  the  averaging  intrinsically,  our  meshless  subdivision  rules  raise  the 
questions  of  how  to  compute  geodesic  centroids  on  manifolds  and  how  to  determine 
a  suitable  neighbour  weighting  scheme.  In  this  paper,  we  guide  the  choice  of  weights 
by  experimental  results  (Section  6)  rather  than  theoretical  evidence  for  the  scheme’s 
convergence  towards  a  smooth  limit  surface.  Future  work  will  consider  formal  proofs 
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Fig.  3.  The  unweighted  centroid  of  a  (blue)  subset  of  this  set  of  points  is  expected  to  be  located 
on  or  near  the  underlying  surface.  Due  to  the  use  of  intrinsic  distances,  this  is  the  case  when 
computing  the  geodesic  centroid  (red).  By  contrast,  in  the  case  of  the  Euclidean  averaging  of 
the  (blue)  points,  the  resulting  centroid  (grey)  is  located  away  from  the  underlying  surface.  This 
effect  gets  more  pronounced  when  increasing  the  size  of  the  subset  (from  left  to  right).  Note  that 
for  geodesically  close  neighbourhoods  and  those  kinds  of  neighbourhoods  only,  the  orthogonally 
projected  {H-m  •  Euclidean  average,  i.e.  the  extrinsic  mean,  generally  provides  a  good 

(first)  estimate  of  the  geodesic  centroid  (leftmost). 


of  the  scheme’s  convergence  to  a  limit  surface  and,  consequently,  any  light  such 
proofs  may  throw  on  the  optimal  choice  of  weights. 

As  regards  the  computation  of  geodesic  centroids,  Buss  and  Fillmore  [2001] 
present  an  algorithm  for  the  computation  of  geodesic  averages  on  spheres.  We 
generalise  the  underlying,  earlier  idea  ([Karcher  1977]  and  references  therein)  of 
minimising  a  least  squares  expression  in  geodesic  distances  in  the  following  section. 

5.3  Geodesic  centroid  computation 

The  benefit  of  performing  the  averaging  intrinsically  is  that  it  ensures  that  sub¬ 
division  generates  smoother,  denser  representations  which  remain  geometrically 
close  to  the  surface.  This  is  not  guaranteed  to  be  the  case  when  considering  Eu¬ 
clidean  instead  of  geodesic  centroids.  For  the  simple  example  illustrated  in  Figure 
3,  Euclidean  averaging  ignores  the  non-linear,  intrinsic  geometry  of  the  object  and 
moves  the  centroid  away  from  the  surface.  By  contrast,  since  the  computation  of 
the  geodesic  centroid  is  based  on  intrinsic  rather  than  Euclidean  distances,  it  is 
inherently  geometry-sensitive  and  falls  onto  the  surface  in  each  case. 

The  weighted  geodesic  centroid  of  a  set  of  n  points  is  defined  as  the  point  g  £  M 
which  minimises  the  weighted  sum  of  squared  intrinsic  distances  to  each  point 

1  ” 

Jig)  ■=  ^^Wk(fM{g,pk), 

^  fc=i 

where  wi,...,Wn  represent  the  point  weights,  with  0  <  <  1, 

In  general,  argmingj(-)  may  not  exist  or  may  not  be  a  single  point.  However,  if 
pi, . . .  ,Pn  are  all  contained  in  a  sufficiently  small  open  geodesic  ball  Bm  on  M, 
a  unique  solution,  pbm  of  J(-))  which  happens  to  lie  in  Bm  [Karcher  1977],  is 
guaranteed.  The  property  we  are  alluding  to  here  is  (geodesic)  convexity,  i.e.  for 
any  Pi^Pj  €  Bm,  the  minimal  geodesic  from  pi  to  Pj  is  unique  in  M  and  contained 
in  Bm- 

In  the  Euclidean  case,  direct  differentiation  of  J(-)  yields  the  minimiser  g  = 
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WkPk-  This  simple  result  does  not  extend  to  the  general  case  considered  here 
but  we  can  prove  that  any  minimiser  must  satisfy: 

n 

yig)  ■=  M^d\i{g,pk)  =  0. 

fc=i 

Then,  starting  from  a  good  initial  guess  g^,  we  can  track  the  minimiser  g  using  back 
propagation  with  velocity  field  V{-).  This  is  due  to  the  fact  that  if  go  G  Bm  and 
Bm  as  above,  then  —V{x)  points  towards  for  x  G  Bm  [Karcher  1977]. 

In  practise,  we  set  go  =  Hm  (X]fe=i  '^kPk),  where  IIm  :  ^  M  is  the  orthog¬ 

onal  projection  operator^.  We  now  show  that  in  the  light  of  the  considerations 
presented  above,  this  extrinsic  mean  represents  either  a  reasonable  initial  condition 
for  the  back  propagation  or  a  first  approximation  to  the  intrinsic  centroid.  By  a 
simple  application  of  Lemma  17  of  [Wallner  and  Dyn  2003],  we  have  that 

/  n 

WkPk  -  IIm  I  X!  '^kPk 

\\k=l  \fc=l 

where  C  is  a  global  constant  which  depends  on  the  curvatures  of  M.  Then,  let 
B  =  BM{x,e),  for  some  x  &  M  and  e  >  0.  Since  \\pi  —  xjj  <  dM{Pi,x)  <  e,  we 
also  have  ||  Ylk=i  ^kPk  -  x\\  <  e.  Therefore,  since  Jj^o  -  a^jj  <  llffo  -  J2k=i  ^kPkW  + 
II  YJk=i  WkPk  -  a:||,  we  obtain 

llffo  ~  a;||  <  +  e, 

which  implies  dMigo^x)  <  e(l  -I-  He){l  +  Ce),  for  another  constant  H  depending 
on  global  metric  properties  of  M  [Memoli  and  Sapiro  2003].  We  only  care  for  a 
simplified  bound 

dnigoix)  <  Et. 

Finally,  let  J  >  0  be  the  maximal  5  >  0  such  that  Bm{x,6)  is  (geodesically) 
convex.  Note  that  it  is  a  fact  that  if  5  <  I  min  {inj{M),  where  inj{M)  is 

the  injectivity  radius  of  M  and  K  bounds  all  sectional  curvatures  in  M  from  above, 
then  Bm{x,6)  is  convex  for  any  x  G  M.  See  §7.6  and  §7.7  in  [Chavel  1997].  For 
such  a  J  >  0  and  provided  e  <  S/E,  and  {pi, . . .  ,p„}  C  Bm{x,  e)  for  some  x  G  M, 
go  G  Bm{x,  5)  and -V (go)  will  be  pointing  towards  gsM-  Also,  in  case  we  want  to  use 
go  as  an  approximation  to  we  have  the  (weak)  bound  dM^gsM  ^  9o)  <  {E+l)t. 
Therefore,  go,  as  defined  above,  represents  a  sensible  choice  as  the  initial  condition 
of  an  eventual  back  propagation  step,  or,  in  any  case,  a  rough  approximation  to 
9Bm  with  known  error  bound.  See  also  Figure  3  (left).  Note  in  particular  that  it  is 
also  a  useful  choice  from  the  point  of  view  of  computational  ease.  The  algorithm  is 
summarised  in  Algorithm  2. 

To  demonstrate  the  applicability  of  this  approach  in  the  context  of  meshless 
subdivision,  we  first  consider  the  case  of  M  representing  the  unit  sphere  in  the 
following  section. 


E 


<  C{diam{B)y 


^To  start  from  any  of  the  points  represents  another  simple  choice. 
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Input:  Intrinsic  Voronoi  neighbourhood  Np  of  point  p  G  Pi.  Weights  Wi  at  points 

Qi  €  Np. 

Output:  Weighted  geodesic  centroid  g. 


0  ***  Computation  of  extrinsic  centroid  go  *** 

1  Compute  the  Euclidean  weighted  centroid  of  Np-, 

2  Compute  go  by  orthogonally  projecting  the  weighted  Euclidean  centroid; 

3 

4  ***  Computation  of  intrinsic  centroid  g  *** 

5  Compute  local  weighted  distance  maps  dn’’  (9i,  •)  from  each  neighbour  qi  G  Np 
outwards  and  accumulate  their  squared  values  at  the  grid  vertices; 

6  Approximate  the  gradient  of  the  accumulated  distance  maps  using  finite 
difference  approximation; 

7  Back  propagate  from  go  towards  g  by  following  the  negative  gradient; 

Alg.  2:  Procedure  for  computation  of  a  weighted  geodesic  centroid  in  pseudocode. 

6.  EXPERIMENTAL  RESULTS  AND  IMPLEMENTATIONAL  DETAILS 

We  begin  with  the  intrinsic  meshless  subdivision  of  a  set  of  points  sampled  relatively 
regularly  uniformly  from  the  surface  of  the  unit  sphere  in  This  initial  restric¬ 
tion  to  spherical  geometry  allows  for  the  computation  of  precise  geodesic  distances 
without  the  need  for  numerical  techniques.  This  way  qualitative  and  quantitative 
aspects  of  our  operator  can  be  presented  without  the  influence  of  the  particular 
projection  and  gradient  descent  techniques  utilised  when  processing  more  complex 
geometry.  This  presentation  is  followed  by  applications  of  our  subdivision  oper¬ 
ator  to  more  complex  geometry.  The  section  concludes  with  comments  on  some 
implementational  aspects. 

6.1  Results  and  discussion 

To  implement  the  intrinsic  centroid  computation  method,  techniques  for  the  com¬ 
putation  of  intrinsic  distances  between  points  on  the  surface,  the  projection  of  the 
starting  point  for  the  back  propagation  onto  the  surface  and  the  computation  of  the 
back  propagation  itself  are  required.  In  the  case  of  the  unit  sphere,  these  techniques 
are  readily  available  and  no  numerical  techniques  are  required.  Intrinsic  distances 
between  points  follow  trigonometrically  and  orthogonal  projection  is  trivial.  Simi¬ 
larly,  the  exponential  map  and  its  inverse  are  directly  available  and  may  be  used  to 
implement  the  back  propagation  procedure.  As  a  result,  for  the  case  of  spherical 
geometry,  our  approach  for  geodesic  centroid  computation  narrows  down  to  the 
technique  of  Buss  and  Fillmore  [2001]. 

We  apply  our  intrinsic  meshless  subdivision  operator  to  a  base  point  set  Pq  of 
2144  points  sampled  relatively  regularly  uniformly  from  the  unit  sphere  (Figure 
5).  The  application  of  our  subdivision  operator  to  Pq  using  the  initial  intrinsic 
proximity  information  from  VD(Po)  yields  the  subdivided  point  set  Pi  of  Figure  6. 
The  result,  P2,  obtained  from  the  application  of  the  operator  to  Pi  using  natural 
neighbourhood  information  from  FD(Pi)  is  shown  in  Figure  7.  Given  the  relatively 
strong  regularity  of  the  input  data,  uniform  weighting  was  used  for  both  the  refine¬ 
ment  and  the  geometric  averaging  rule  in  both  iterations.  The  results  produced  by 
our  meshless  subdivision  operator  are  presented  alongside  the  point  sets  produced 
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Fig.  4.  Histograms  of  the  (spherical)  distance  from  each  point  in  Pi  (left)  and  P2  (right)  to  its 
closest  neighbour  in  the  respective  set. 


Model\k 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Pi 

0.807 

0.832 

0.832 

0.821 

0.788 

0.809 

0.8132 

0.786 

0.818 

0.828 

P2 

0.781 

0.815 

0.823 

0.804 

0.785 

0.799 

0.7980 

0.784 

0.810 

0.822 

Table  I.  Values  of  the  density  uniformity  measure  p{k)  for  Pi  and  P2  and  with  k  £  {1,  2, . . . ,  10}. 
The  values  underline  the  regular  uniformity  of  the  subdivided  point  sets  generated  by  our  algo¬ 
rithm. 

by  the  application  of  Loop  subdivision  to  a  triangular  mesh  representation  of  Pq- 
As  indicated  by  the  detail  views  of  Figure  7,  the  point  distributions  obtained  from 
the  two  operators  after  two  iterations  are  qualitatively  similar  with  the  slight  ir¬ 
regularities  in  the  distribution  of  P2  being  slightly  more  pronounced  in  the  case  of 
meshless  subdivision  due  to  the  use  of  uniform  weights.  There  are,  however,  no 
noticeable  differences  in  the  smoothing  effect  of  these  two  operators. 

In  order  to  analyse  the  point  set  distributions  generated  by  our  subdivision  opera¬ 
tor  quantitatively,  we  compute  the  mean  and  the  standard  deviation  of  the  distance 
from  each  point  in  the  set  to  its  closest  neighbour  (s)  for  the  subdivided  point  sets  Pi 
and  P2-  For  each  p  in  the  point  set  Pi,  let  sdk{p)  denote  the  (spherical)  distance  from 
p  to  its  kth  closest  neighbour.  As  an  indicator  for  the  uniformity  of  the  density  of 
point  set  Pi,  consider  p{k)  =  sdt (p)  •  Since  p{k)  represents  an  absolute  measure, 

it  may  be  too  sensitive,  therefore  we  compute  instead  p{k)  =  ™gan(sd'‘)+std(sd'')  ’ 
where  mean  and  std  stand  for  the  mean  and  standard  deviation  of  the  spherical 
distances  over  the  point  set  respectively. 

The  histograms  of  sdi{x)  corresponding  to  the  two  sets  of  points  are  given  in 
Figure  4.  Note  in  particular  in  Table  I  that  the  values  of  p{k),  for  1  <  A:  <  10, 
are  quite  close  to  1.0  therefore  indicating  small  dispersion  up  to  the  10th  closest 
neighbour. 

Figure  8  presents  an  application  example  dealing  with  more  complex  geometry. 
We  apply  our  meshless  subdivision  operator  to  a  base  point  set  of  10088  points 
generated  from  the  Michelangelo  Youthful  data  set  with  the  help  of  [Moenning  and 
Dodgson  2004] .  Experimentation  revealed  the  base  point  set  to  be  regular  enough  to 
allow  for  simple  reciprocal  distance  weighting  in  the  computations  of  the  weighted 
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geodesic  centroids.  The  flatly  shaded  renderings  of  the  surfaces  reconstructed  from 
the  subdivided  point  sets  Pi  and  P2  illustrate  the  smoothing  effect  of  the  meshless 
subdivision.  As  indicated  by  the  comparative  illustrations  in  Figure  9,  this  approach 
may  be  used  to  obtain  a  smoother,  more  regular  and  more  compact  representation 
of  an  highly  dense  point  cloud.  A  similar  effect  is  shown  in  Figure  10.  The  50% 
decimated  versions  of  the  rocker  arm  and  screwdriver  CAD  data  sets  available 
from  the  Cyberware  website  were  meshlessly  subdivided  twice.  The  smoothing 
effect  of  these  iterations  is  again  clearly  apparent  when  comparing  the  surfaces 
reconstructed  from  the  subdivided  point  sets  to  those  reconstructed  from  Pq  and 
the  non-simplifled,  non-subdivided  data  sets  respectively. 

The  detail  view  of  Figure  8  highlights  the  local  clustering  effect  caused  by  over¬ 
lapping  neighbour  relations  as  discussed  in  Section  5.1,  an  effect  typically  only 
encountered  with  point  neighbourhoods  which  determine  the  neighbour  relations 
independently  of  each  other.  Our  discrete  approximation  of  the  intrinsic  Voronoi 
diagram,  however,  implies  discretisation  error  and  thus  individual  instances  of  a 
point  being  assigned  to  the  wrong  Voronoi  region.  The  generally  limited  number 
of  these  instances  is  considered  preferable  to  the  complications  associated  with  ad¬ 
dressing  the  problem  of  overlapping  neighbour  relations  (see,  e.g.  Guennebaud  et  al. 
[2004])  when  using  alternative  neighbourhood  concepts. 

The  limited  impact  of  these  assignment  errors  is  illustrated  by  the  detail  view 
of  Figure  11.  The  Isis  data  set  was  subdivided  once  with  the  regularity  of  the 
subdivided  point  set  being  only  mildly  affected  by  erroneous  neighbourhood  as¬ 
signments.  In  contrast  to  the  processing  of  the  Youthful  base  point  set,  meshless 
subdivision  of  the  Isis  point  cloud  using  geodesic  centroids  was  found  to  yield  results 
not  significantly  different  from  the  more  efficient  subdivision  by  extrinsic  centroids. 
We  exploit  this  observation,  due  to  the  high  density  of  the  initial  point  set,  by 
presenting  the  results  from  orthogonal  projection  of  the  uniformly  weighted  Eu¬ 
clidean  centroids,  i.e.  without  subsequent  gradient  descent  towards  their  geodesic 
counterparts. 

Table  II  summarises  the  parameter  settings  and  point  set  sizes  for  the  various 
application  examples.  Using  our  non-optimised  implementation,  the  offset  band 
computation  pre-processing  step  and  a  meshless  subdivision  iteration  took  a  maxi¬ 
mum  of  a  few  hundred  seconds  each.  Both  offset  band  and  intrinsic  Voronoi  diagram 
computation  efficiency  generally  depend  strongly  on  the  offset  band  radius  r  and  the 
grid  spacing,  our  settings  of  which  are  presented  in  the  table.  Note  in  this  context 
that  the  numerical  error  inherent  in  the  intrinsic  distance  computations  increases 
with  r,  the  admissible  minimum  value  of  which  necessarily  increases  with  the  grid 
spacing  [Memoli  and  Sapiro  2001].  Details  on  various  aspects  of  the  implementation 
of  our  algorithm  are  provided  in  the  following  section. 

6.2  Implementational  details 

The  algorithm  was  implemented  in  C-I--I-  (Microsoft  Visual  C-|— I-  7.1)  with  the  help 
of  the  “Blitz-|— I-  0.7  Numerical  Library”  [Blitz-|--|-  2003]  on  a  Pentium  4  2.8GHz, 
512MB  Windows  XP  machine.  In  the  following,  we  summarise  relevant  implemen¬ 
tational  details. 
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Iteration  1 


Ax(=  Ay  =  Az) 

r 

Xol 

\Pi\ 

Sphere 

0.25 

0.8 

2144 

8570 

Youthful 

1.0 

2.2 

10080 

39888 

Screwdriver 

0.5 

1.0 

13577 

56220 

Rocker  arm 

1.0 

1.9 

20088 

81442 

Isis 

0.1 

0.2 

187644 

760162 

Iteration  2 


Ax(=  Ay  =  Az) 

r 

\Pi\ 

\P2\ 

Sphere 

0.25 

0.75 

8570 

34275 

Youthful 

0.25 

0.75 

39888 

208010 

Screwdriver 

0.25 

0.6 

56220 

295110 

Rocker  arm 

0.25 

0.45 

81442 

488212 

Table  II.  Parameter  settings  and  point  set  sizes  for  the  application  examples;  Ax,  Ay,  Az  refer  to 
the  grid  spacing  in  the  three  principal  Cartesian  grid  directions;  r  represents  the  constant  offset 
band  radius. 

6.2.1  Offset  band  computation  and  intrinsic  distance  mapping.  Intrinsic  dis¬ 
tance  mapping  requires,  firstly,  the  computation  of  the  offset  band  flp  and,  sec¬ 
ondly,  weighted  Euclidean  distance  mapping  within  the  offset  band.  We  meet  both 
of  these  requirements  with  the  help  of  conventional  Fast  Marching.  For  a  given 
r,  offset  band  computation  amounts  to  the  simultaneous  propagation  of  (circular) 
fronts  from  each  input  point  outwards  until  the  front’s  extent  equals  r.  The  off¬ 
set  band  consists  of  those  grid  vertices  visited  during  the  propagation.  All  other 
grid  vertices  of  the  point  set’s  bounding  volume  are  discarded  to  minimise  memory 
usage.  When  dealing  with  relatively  large  point  sets,  we  perform  this  offset  band 
computation  as  a  pre-processing  step  which  makes  the  set  of  valid  grid  vertices 
available  for  subsequent  weighted  distance  mapping.  Irregularly  distributed  point 
sets  may  require  the  offset  band  to  be  adaptive.  This  can  be  achieved  by  computing 
radii  using  a  minimum  spanning  tree  or  local  principal  component  analysis  as 
discussed  in  [Memoli  and  Sapiro  2001;  2003].  Weighted  distance  mapping  within 
flp  (or  flp)  is  achieved  by  another  application  of  conventional  Fast  Marching  re¬ 
stricted  to  operate  within  the  offset  band  only,  i.e.  the  set  of  grid  vertices  returned 
by  the  preceding  band  fitting  step. 

A  single  min-heap  as  typically  used  for  the  efficient  implementation  of  the  con¬ 
ventional  Fast  Marching  method  [Sethian  1999]  and  a  Cartesian  grid  represent  the 
main  data  structures  required  in  this  context.  The  grid  data  structure  is  imple¬ 
mented  in  the  form  of  a  lookup  table  with  each  grid  vertex  mapped  to,  amongst 
others,  its  distance  map  value  and  min-heap  index  and  offset  band  membership 
status.  By  using  a  lookup  table,  invalid  grid  vertices  can  be  quickly  discarded  and 
different  grid  spacing  in  each  direction  is  supported. 

6.2.2  Geodesic  Voronoi  diagrams  and  intrinsic  natural  neighbourhoods.  The  nat¬ 
ural  neighbourhood  information  of  a  point  Pi  G  Pi  is  available  from  VD(Pi).  We 
approximate  this  diagram  using  weighted  geodesic  distance  maps,  i.e.  in  analogy  to 
the  dropping  of  pebbles  in  still  water,  circular  fronts  move  across  the  surface  from 
the  points  of  impact.  The  locations  where  wave  fronts  meet  define  the  geodesic 
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Voronoi  diagram  of  the  points  of  impact.  Figure  12  gives  a  triangular  mesh-based 
example  produced  using  public  domain  software  [Peyre  and  Cohen  2003].  The 
wave  propagation  is  discretised  and  simulated  accurately  by  solving  (2)  using  our 
implementation  of  the  intrinsic  distance  mapping  for  surfaces  in  point  cloud  form 
discussed  in  Section  6.2.1.  This  way  Voronoi  edges/ vertices  and  thus  each  point’s 
Voronoi  neighbours  are  obtained  during  front  propagation  as  loci  of  intersection 
between  geodesic  offset  curves  [Cohen  2001]. 

The  initial  geodesic  Voronoi  diagram,  VD(Po),  is  either  computed  during  the 
algorithm’s  initialisation  phase  by  propagating  fronts  simultaneously  from  all  pi  G 
Pq  outwards  (Figure  12)  or  follows  from  prior  intrinsic  point  cloud  simplification  as 
indicated  in  Section  5.2. 

The  neighbourhood  information  is  held  in  the  form  of  a  lookup  table  mapping 
each  p  to  its  set  of  natural  neighbours  Np  represented  by  a  set  of  indices  referencing 
the  corresponding  input  points.  The  lookup  table  of  (valid)  grid  vertices,  already 
discussed  above,  represents  the  only  other  main  data  structure  used  in  this  context 
and  is  required  here,  as  above,  for  distance  mapping  support. 

6.2.3  Geodesic  centroid  computation.  Using  the  natural  neighbourhood  infor¬ 
mation  Np  of  point  p,  its  weighted  Euclidean  centroid  is  readily  available.  We 
use  the  “almost”  orthogonal  projection  operator  of  Alexa  and  Adamson  [2004]  to 
project  this  centroid  onto  M  thereby  obtaining  g^  (Section  5.3).  Fast  Marching 
for  point  clouds  [Memoli  and  Sapiro  2003]  is  then  again  used  to  compute  distance 
maps  within  the  offset  band  from  each  point  in  Np  outwards  with  the  squared  dis¬ 
tance  map  values  being  accumulated  at  the  grid  vertices.  To  avoid  unnecessary 
propagation,  the  extent  of  the  distance  mapping  is  limited  to  the  radius  of  the 
sphere  enclosing  the  points  in  Np.  A  standard  Runge-Kutta  gradient  descent  pro¬ 
cedure  with  multilinear  interpolation  is  finally  employed  to  back  propagate  from 
(/o  towards  the  weighted  geodesic  centroid  g  by  following  the  (negative)  gradient  of 
the  accumulated  distance  maps  estimated  with  the  help  of  the  normalised  central 
difference  operator  provided  by  Blitz-|— I-  0.7  [Blitz-|— I-  2003]. 

7.  CONCLUSION 

We  introduced  the  concept  of  meshless,  or  point  cloud,  subdivision  based  on  the 
computation  of  weighted  geodesic  centroids  on  manifolds  represented  by  noise-free 
or  noisy  point  clouds.  We  implemented  and  showed  the  applicability  of  this  tech¬ 
nique  with  the  help  of  a  new  method  for  the  computation  of  such  centroids. 

The  consideration  of  local  intrinsic  point  proximity  instead  of  mesh  connectivity 
helps  to  overcome  some  of  the  limitations  of  existing  mesh  subdivision  schemes. 
For  example,  by  working  with  the  raw  data  and  operating  directly  across  the  point 
cloud,  problems  associated  with  the  complexity  of  the  topological  subdivision  of 
high-dimensional  meshes  are  avoided.  For  extremely  high-dimensional  data,  mesh¬ 
less  subdivision  operators  need  to  be  devised,  however,  which  upsample  the  point- 
sampled  geometry  more  conservatively  than  the  operator  suggested  in  this  paper. 
In  this  direction,  we  are  considering  introducing  adaptive  neighbourhoods  based  on 
curvature  estimators  such  as  those  reported  in  Cazals  and  Pouget  [2003]  and  Mitra 
and  Nguyen  [2003].  We  leave  this  to  future  work. 

As  first  proposed  in  the  context  of  univariate  spline  subdivision  schemes  [Lane 
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and  Riesenfeld  1980],  we  plan  to  investigate  the  practical  and  theoretical  aspects 
of  repeating  the  geometric  averaging  step  several  times  after  each  refinement  step. 
We  expect  to  get  higher  smoothness  with  a  higher  number  of  averaging  steps  in 
each  iteration  of  the  meshless  subdivision  process. 

Depending  on  the  extent  of  any  non-uniformities  of  the  input  point  cloud,  the 
experimental  selection  of  point  weights  tends  to  be  elaborate.  We  are  currently 
working  on  the  more  systematic  choice  of  weights.  This  issue  is  of  course  closely 
related  to  the  theoretical  analysis  of  our  meshless  intrinsic  subdivision  scheme. 
In  this  context,  we  intend  to  build  upon  ideas  of  Wallner  and  Dyn  [Wallner  and 
Dyn  2003;  2004;  2005],  which  this  research  was  inspired  by,  in  particular  their 
convergence  and  smoothness  analysis  of  non-linear  geodesic  curve  subdivision  by 
proximity  to  a  corresponding  linear  extrinsic  subdivision  scheme. 
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Fig.  5.  Relatively  regularly  uniformly  distributed  base  point  set,  Pq,  of  2144  points 
acquired  from  the  unit  sphere  (right).  The  triangular  base  mesh  generated 
from  Pq  for  the  support  of  Loop  subdivision  is  shown  on  the  left.  A  flatly 
shaded  rendering  of  the  reconstructed  surface  is  shown  in  the  centre. 


Fig.  6.  Results  after  one  iteration  of  Loop  (left)  and  meshless  subdivision  (right); 
|Pi|  =  8570  points  in  both  cases.  Flatly  shaded  renderings  of  the  recon¬ 
structed  surfaces  are  shown  next  to  each  point  set. 


Fig.  7.  Results  after  the  second  iteration  of  Loop  (left;  34275  points)  and  mesh¬ 
less  subdivision  (right;  36442  points).  The  corresponding  reconstructed 
surfaces  are  given  next  to  each  point  set.  The  detail  views  indicate  how 
the  slight  irregularities  in  the  distribution  of  Pq  lead  to  slightly  more  pro¬ 
nounced  local  irregularities  in  the  case  of  meshless  subdivision  due  to  the 
use  of  uniform  weighting.  These  slightly  more  pronounced  irregularities  do 
not  have  any  noticeable  effect  on  the  smoothness  of  the  surface. 
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Fig.  8.  Flatly  shaded  renderings  of  the  surfaces  reconstructed  from  point  sets  l^ol  = 
10088  (top  and  bottom  left),  |Pi|  =  39888  (centre  and  bottom  right)  and 
IP2I  =  208010.  The  smoothing  effect  of  meshless  subdivision  is  clearly  visible. 
As  illustrated  by  the  detail  front  and  side  views,  the  subdivided  point  sets 
are  distributed  relatively  regularly  uniformly.  Instances  of  local  irregularities 
(encircled  in  black)  are  caused  by  discretisation  errors  during  discrete  intrinsic 
Voronoi  diagram  computation. 
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Fig.  9.  Flatly  shaded  renderings  of  the  surfaces  reconstructed  from  the  non-simplified 
Youthful  point  set  (left)  and  P2  of  Figure  8  illustrating  how  meshless  subdivi¬ 
sion  of  a  base  point  set  resulting  from  point  cloud  simplification  may  be  used 
to  obtain  a  smoother,  more  regular  and  more  compact  representation  of  the 
original  data  set  of  1728305  points. 


Fig.  10.  On  the  left,  the  flatly  shaded  renderings  of  surfaces  reconstructed  from  the 
50%  decimated  screwdriver  (top)  and  rocker  arm  (bottom)  CAD  data  sets 
(Pq)  are  shown.  The  reconstructed  surfaces  obtained  from  two  iterations  of 
meshless  subdivision  of  these  point  sets  are  presented  in  the  centre  (^2)-  The 
surfaces  reconstructed  from  the  non-simplified,  non-subdivided  CAD  data  sets 
are  given  on  the  right.  The  smoothing  effect  of  meshless  subdivision  is  again 
clearly  visible. 
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Fig.  11.  Point  set  detail  views  and  smoothly  shaded  renderings  of  the  Isis  data  set  (left), 
Po)  I-PqI  =  187644,  and  the  subdivided  point  set  Pi  (right),  |Pi|  =  760162, 
resulting  from  one  iteration  of  meshless  subdivision.  Due  to  the  high  density  of 
Pq,  the  difference  in  location  between  the  (uniformly  weighted)  extrinsic  and 
the  corresponding  intrinsic  centroid  was  found  to  be  negligible  and  extrinsic 
means  were  used  throughout. 


Fig.  12.  Wave  propagation  for  the  computation  of  a  discrete  geodesic  Voronoi  diagram. 
By  simultaneously  propagating  waves  for  geodesic  distance  mapping  purposes 
from  the  (red)  generator  points  outwards  (left),  an  intrinsic  Voronoi  partition¬ 
ing  of  the  triangulated  surface  is  obtained  (right). 
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